Librerías

library(geoR)
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------
library(geoR,lib.loc="PATH_TO_geoR")
library(sp)
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tmap)
#> Error : invalid version specification '1,5'
library(e1071)
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(gstat)
library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
library(nlme)
library(ade4)
## 
## Attaching package: 'ade4'
## The following object is masked from 'package:spdep':
## 
##     mstree

Contexto

Este conjunto de datos contiene el contenido de magnesio en la capa de suelo de 20 a 40 cm en 178 ubicaciones dentro de una determinada área de estudio dividida en tres subáreas. La elevación en cada ubicación también se registró.La primera región suele estar inundada durante la temporada de lluvias y no se usa como área experimental. Los niveles de calcio representarían el contenido natural en la región. La segunda región ha recibido fertilizantes hace un tiempo y está típicamente ocupada por campos de arroz. La tercera región ha recibido fertilizantes recientemente y se usa con frecuencia como área experimental.

Formato

El set de datos original es un data frame con 178 observaciones y 10 variables, de las cuales para este trabajo solo son consideradas las siguientes:

Estos datos están disponibles en el paquete geoR (Paulo J. Ribeiro Jr and Peter J. Diggle (2016)).

Para realizar los análisis exploratorios a seguir, utilizaremos el paquete geoR e inicialmente necesitaremos crear un objeto del tipo geodata, como sigue.

Recoleción de datos

Importamos los datos para R

datos <- read.csv("~/Regresion Avanzada/Grupo10.csv", header=TRUE)
View(datos)
# Mostrar de forma compacta la estructura de un objeto R arbitrario
str(datos)  
## 'data.frame':    178 obs. of  4 variables:
##  $ X     : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ east  : int  5710 5727 5745 5764 5781 5799 5817 5837 5873 5907 ...
##  $ north : int  4829 4875 4922 4969 5015 5062 5109 5161 5254 5347 ...
##  $ mg2040: int  16 21 22 12 36 22 25 25 28 18 ...
head(datos)
##   X east north mg2040
## 1 1 5710  4829     16
## 2 2 5727  4875     21
## 3 3 5745  4922     22
## 4 4 5764  4969     12
## 5 5 5781  5015     36
## 6 6 5799  5062     22
class(datos)
## [1] "data.frame"

Transformamos el data set a un objeto de geodata

datos_geo <- as.geodata(datos, coords.col=2:3, data.col=4)
class(datos_geo)
## [1] "geodata"

Para hacer el análisis de datos espaciales de manera más práctica, es común transformar el sistema de coordenadas geográficas en un sistema de coordenadas cartesianas, como el sistema de proyección UTM. Esto permite que las distancias entre los sitios se expresen en metros en lugar de grados, lo que facilita el análisis. Por lo tanto, el primer paso en el análisis de datos espaciales sería convertir las coordenadas geográficas en coordenadas cartesianas (UTM), en este caso no es necesario dado que las variables north y east que representan las coordenadas se encuentran expresadas en metros.

Análisis Exploratorio Espacial

En el análisis exploratorio de datos geoestadísticos es importante explorar la distribución de la variable. Esto se puede hacer mediante estadísticas descriptivas y gráficos de distribución de frecuencias. Si se supone que los datos tienen una distribución normal, estas medidas pueden ayudar a verificar los supuestos. Una distribución simétrica y similar a la de una variable normal se da cuando la media y la mediana son casi iguales y el coeficiente de asimetría es inferior a 1. La distribución de la variable también puede ayudar en la depuración de datos raros.

Breve resumen descriptivo de los datos

#Da un rápido resumen del objeto geodata cargado.
summary(datos_geo)
## Number of data points: 178 
## 
## Coordinates summary
##     east north
## min 4957  4829
## max 5961  5720
## 
## Distance summary
##        min        max 
##   43.01163 1138.11774 
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 10.00000 21.00000 25.50000 25.58427 30.00000 43.00000

Gráficamos los valores en las posiciones espaciales (distinguiendo según cuartiles), los datos frente a las coordenadas y un histograma de los datos. Los gráficos de dispersión de los datos frente a las coordenadas nos pueden ayudar a determinar si hay una tendencia

plot(datos_geo)

Generamos un gráfico con las posiciones de los datos (y por defecto con el tamaño de los puntos proporcional al valor):

 points(datos_geo, col = "gray", pt.divide = "equal")

Transformación

datos_sf <- st_as_sf(datos, coords=c("east","north"), 
                     crs=32721)
head(datos_sf)
## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 5710 ymin: 4829 xmax: 5799 ymax: 5062
## Projected CRS: WGS 84 / UTM zone 21S
##   X mg2040          geometry
## 1 1     16 POINT (5710 4829)
## 2 2     21 POINT (5727 4875)
## 3 3     22 POINT (5745 4922)
## 4 4     12 POINT (5764 4969)
## 5 5     36 POINT (5781 5015)
## 6 6     22 POINT (5799 5062)

Eliminación de outliers e inliers

Las funciones hist() y boxplot() realizan gráficos de histogramas y box-plots, respectivamente. Sus múltiples argumentos permiten la edición de cada gráfico

par(mfrow = c(1, 2))
hist(
  datos_sf$mg2040,
  col = 'grey',
  nclass = 20,
  main = "Histograma",
  ylab = 'Frecuencia Relativa',
  xlab = 'Contenido Magnesio (mmolc/dm3)'
)

boxplot(
  datos_sf$mg2040,
  col = 'grey',
  ylab = 'Contenido Magnesio (mmolc/dm3)',
  main = "Box-Plot",
  ylim = c(0, 50)
)

La función skewness() del paquete e1071 permite calcular el coeficiente de asimetría. Existen 3 fórmulas para su cálculo (por defecto usa el tipo 3)

skewness(datos_sf$mg2040)
## [1] 0.1099413

En el histograma se observa que los datos tienen una distribución similar a una normal. No se ve una asimetría hacia algun lado en particular sino una disribución uniforme, también puede advertirse con los estadísticos con el coeficiente de asimetría el cual es de 0,109. En el gráfico box-plot no se observan valores extremos de la variable que se encuentren principalmente por encima de la media + 3SD

Las siguientes instrucciones calculan y crean objetos para la media, el DE y los límites superior (media + 3DE) e inferior (media - 3DE) con los que pueden detectarse los outliers y eliminar restos valores caso que sea necesario.

Media <- mean(datos_sf$mg2040)
DE <- sd(datos_sf$mg2040)
LI <- Media-3*DE
LS <- Media+3*DE
# Rango de valores de outliers
(datos_sf$mg2040 > LS | 
           datos_sf$mg2040 < LI)
##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Como no existe ningun valor dentro de este rango no es necesario la eliminación de outlieres

Detección de Inliers

La identificación y eliminación de inliers requiere de la creación de una matriz de ponderación espacial. La función dnearneigh() se utiliza para identificar el vecindario de cada sitio. Para ello es necesario calcular la distancia espacial entre los puntos para lo cual se usa la sintaxis $geom que permite acceder a las coordenadas del objeto datos_sf. En este ejemplo, se consideran datos vecinos a aquellos que se encuentran a una distancia Euclídea de 0 a 100 m. La función nb2listw() transforma el objeto vecindarios que contiene las distancias a una matriz de pesos estandarizados por filas (style = "W"). Este objeto es denominado pesos_sp. Para generar la matriz de pesos espaciales es necesario que todos los puntos tengan al menos un vecino, caso contario la función nb2listw() genera un error advirtiendo este hecho.

# 8345 es de hacer la distancia norte y este al cuadrado raiz seria el valor máximo
vecindarios <- dnearneigh(datos_sf$geom,
                          d1 = 0, d2 = 100)

# Muestra la cantidad de vecinos para la posicion
card(vecindarios)
##   [1]  3  5  7  6  6  7  6  6  5  6  7  5  5  6  8  8  9  9  9  7  7  7  8  8  9
##  [26]  6  7  7  9 10 10 11 10  8  8  9  8  7  9 10 11 10 10 11  9  9  9  8  8  5
##  [51]  6  7 10 11  9 10  9  8  8  8  8  5  5  7 10 12 11 10 10 10  9  7  7  6  7
##  [76] 10 12 11 10  9  9  8  7  5  3  3  8 10 12 11 10 10 10 10  8  4  2  9 10 12
## [101] 11 11 10 10  9  9  7  5  6  8 10 10 10  9  9  6  7  5  8 11 11 11 10  9 10
## [126]  6  4  2  5 10 10 11 10 10 10 10  6  3  4  7 10 11 10 12 10  9  9  4  3  5
## [151]  5  6 10 11 11 10  9  9  8  6  1  5  2  7 10 10 10  8  7  5  1  3  3  6  8
## [176]  8  7  6
# Hace un resumen de todos los valores
summary(vecindarios)
## Neighbour list object:
## Number of regions: 178 
## Number of nonzero links: 1416 
## Percentage nonzero weights: 4.469133 
## Average number of links: 7.955056 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 12 
##  2  3  7  4 15 17 21 22 25 41 16  5 
## 2 least connected regions:
## 161 171 with 1 link
## 5 most connected regions:
## 66 77 89 100 144 with 12 links
pesos_sp <- nb2listw(vecindarios,
                     style = "W")

La función localmoran() calcula el índice local de Moran que permite identificar potenciales outliers espaciales. La información referida al valor del índice local de Moran de cada punto se encuentra en la columna Ii mientras que su significancia estadística en la columna Pr(z < 0).

moranl <-
  localmoran(datos_sf$mg2040,
             pesos_sp,
             alternative = "less")
head(moranl)
##            Ii         E.Ii     Var.Ii       Z.Ii Pr(z < E(Ii))
## 1  0.86581506 -0.011967143 0.69358099  1.0539941     0.8540572
## 2  0.48460289 -0.002737869 0.09499218  1.5812073     0.9430847
## 3  0.39070018 -0.001673685 0.04103973  1.9368582     0.9736187
## 4 -1.06995511 -0.024040584 0.67628376 -1.2718378     0.1017154
## 5 -0.42053851 -0.014133570 0.40162630 -0.6412803     0.2606703
## 6 -0.08158934 -0.001673685 0.04103973 -0.3944842     0.3466118

La identificación y eliminación de valores atípicos requiere la creación de una matriz de ponderación espacial. Para ello, se utiliza la función “dnearneigh()” para identificar el vecindario de cada sitio. Esto se logra calculando la distancia espacial entre los puntos, para lo cual se utiliza la sintaxis “$geom”, que permite acceder a las coordenadas del objeto “datos_sf”.

En este ejemplo, se consideran datos vecinos a aquellos que se encuentran a una distancia Euclidiana de entre 0 y 100 metros. La función “nb2listw()” transforma el objeto “vecindarios”, que contiene las distancias, en una matriz de pesos estandarizados por filas (style = “W”). Este objeto se denomina “pesos_sp”.

Es importante tener en cuenta que, para generar la matriz de pesos espaciales, es necesario que todos los puntos tengan al menos un vecino. En caso contrario, la función “nb2listw()” generará un error para advertir sobre este hecho.

moranp <-
  moran.plot(
    datos_sf$mg2040,
    col = 3,
    pesos_sp,
    labels = F,
    quiet = T,
    xlab = "Contenido Magnesio",
    ylab = "Contenido Magnesio Spatially Lagged"
  )

Para visualizar en una tabla los puntos potencialmente influyentes y sus estadísticos de diagnóstico se imprime el objeto moranp. Datos con * en la columna inf se los considera como influyente y por lo tanto posible outlier espacial.

summary(moranp)
##        x               wx          is_inf           labels         
##  Min.   :10.00   Min.   :10.00   Mode :logical   Length:178        
##  1st Qu.:21.00   1st Qu.:23.38   FALSE:167       Class :character  
##  Median :25.50   Median :25.65   TRUE :11        Mode  :character  
##  Mean   :25.58   Mean   :25.77                                     
##  3rd Qu.:30.00   3rd Qu.:28.50                                     
##  Max.   :43.00   Max.   :32.30                                     
##      dfb.1_               dfb.x                dffit                cov.r      
##  Min.   :-0.6807993   Min.   :-0.3462534   Min.   :-0.6992421   Min.   :0.852  
##  1st Qu.:-0.0229104   1st Qu.:-0.0242642   1st Qu.:-0.0474738   1st Qu.:1.006  
##  Median :-0.0022348   Median : 0.0011303   Median :-0.0012518   Median :1.016  
##  Mean   :-0.0001416   Mean   : 0.0001139   Mean   : 0.0001092   Mean   :1.011  
##  3rd Qu.: 0.0247714   3rd Qu.: 0.0242995   3rd Qu.: 0.0733676   3rd Qu.:1.019  
##  Max.   : 0.3771651   Max.   : 0.6195409   Max.   : 0.3847945   Max.   :1.059  
##      cook.d               hat          
##  Min.   :5.000e-08   Min.   :0.005640  
##  1st Qu.:3.711e-04   1st Qu.:0.006374  
##  Median :2.449e-03   Median :0.008144  
##  Mean   :6.454e-03   Mean   :0.011236  
##  3rd Qu.:6.341e-03   3rd Qu.:0.013070  
##  Max.   :2.227e-01   Max.   :0.044910

Desde el objeto moranp puede extraerse una matriz de valores lógicos (verdadero/falso) para los estadísticos diagnóstico que identifican un punto como influyente o no.

influ <- moranp$is_inf
head(influ)
## [1] FALSE FALSE FALSE  TRUE FALSE FALSE

En la siguiente sentencia se adiciona al objeto datos_sf los valores de los objetos moranl e influ, que tienen información para detectar los outliers espaciales detectados con el índice de Moran local y gráfico de Moran, respectivamente.

datos_sf_1 <- cbind(datos_sf, moranl, influ)

Posteriormente procedemos a eliminar los datos con Índice de Moran Local negativo y estadísticamente significativos (p<0,05). La función subset() selecciona datos que cumplen con cierta condición lógica. El operador lógico or que indica que extraiga los datos que cumplen con alguna de las dos condiciones. El nuevo objeto es denominado datos_sf_2. Además, se crea un nuevo objeto que tendrá los datos que han sido eliminados en este proceso (inliers_ml).

datos_sf_2 <-
  subset(datos_sf_1, 
         datos_sf_1[["Ii"]] >= 0 | 
           datos_sf_1[["Pr.z...E.Ii.."]] > 0.05)
inliers_ml <-
  subset(datos_sf_1, 
         datos_sf_1[["Ii"]] < 0 &
           datos_sf_1[["Pr.z...E.Ii.."]] < 0.05)
summary(datos_sf_2$mg2040)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.00   21.00   25.50   25.56   30.00   43.00
skewness(datos_sf_2$mg2040)
## [1] 0.1121928
par(mfrow = c(1, 2))
hist(
  datos_sf_2$mg2040,
  col = 'grey',
  nclass = 20,
  main = "Histograma",
  ylab = 'Frecuencia Relativa',
  xlab = 'Contenido Magnesio (mmolc/dm3)'
)
boxplot(
  datos_sf_2$mg2040,
  col = 'grey',
  ylab = 'Contenido Magnesio (mmolc/dm3)',
  main = "Box-Plot",
  ylim = c(0, 50)
)

Luego de identificar y eliminar los inliers detectados con el índice de Moran y posteriormente con el gráfico de Moran, la nueva base de datos presenta 166 casos, es decir, se eliminaron 12 casos (6.7% de los datos) respecto a la base sin outliers.

Los estadísticos descriptivos de los datos depurados muestran una mejora en el coeficiente de asimetría (-0,19) lo cual se refleja en el histograma y box-plot.

Índice de Moran

Para la conformación de la matriz de ponderadores espaciales se definieron los vecindarios de cada sitio mediante una red de conexión construida en base a la distancia euclídea. Se consideraron sitios vecinos a aquellos contiguos ubicados hasta 100 m de distancia. El procedimiento es similar al empleado para el cálculo de índice de Moran local. En este caso se agrega el argumento zero.policy=T dentro de la función nb2listw() y moran.mc(). Esto permite que se genera la matriz de pesos espaciales sin la restricción de que todos los puntos tengan al menos un dato vecino.

vecindarios <- dnearneigh(datos_sf_2$geom, 
                          0, 100)
pesos_sp <- nb2listw(vecindarios, 
                     style = "W", 
                     zero.policy = TRUE)

Para realizar el cálculo del Índice de Moran y determinar su significancia estadística mediante simulación Monte Carlo, se utiliza moran.mc(). Se debe especificar la variable en estudio, la lista con los pesos de las ponderaciones espaciales y el número de simulaciones.

i.moran <-
  moran.mc(
    datos_sf_2$mg2040,
    listw = pesos_sp,
    nsim = 999,
    zero.policy = T
  )
i.moran
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  datos_sf_2$mg2040 
## weights: pesos_sp  
## number of simulations + 1: 1000 
## 
## statistic = 0.37802, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Estos resultados permiten concluir que existe autocorrelación espacial positiva (0,37802) y que esta es estadísticamente significativa (p=0,001).

Índice de Geary

geary.test(datos_sf_2$mg2040, listw = pesos_sp)
## 
##  Geary C test under randomisation
## 
## data:  datos_sf_2$mg2040 
## weights: pesos_sp 
## 
## Geary C statistic standard deviate = 8.6896, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.612803101       1.000000000       0.001985465

El índice de Geary C es una medida de autocorrelación espacial que se utiliza comúnmente en análisis espaciales y de patrones espaciales. El resultado muestra que el valor de Geary C calculado para los datos es significativamente menor que el valor esperado (1.000000000), lo que indica una fuerte autocorrelación espacial en los datos. El valor de Geary C es 0.612803101, lo que sugiere que los valores similares están agrupados juntos en el espacio, mientras que los valores disímiles están más dispersos.

El estadístico desviación estándar de Geary C (8.6896) es alto y el valor p es extremadamente bajo (<2.2e-16), lo que indica que la probabilidad de obtener un valor tan bajo de Geary C por azar es extremadamente baja. Por lo tanto, se rechaza la hipótesis nula de que no hay autocorrelación espacial en los datos y se acepta la hipótesis alternativa de que existe una autocorrelación espacial significativa en los datos.

En resumen, el resultado del índice de Geary C sugiere que los valores similares están agrupados juntos en el espacio y que hay una fuerte autocorrelación espacial en los datos.

Análisis basado en semivariogramas opcion 1

Las semivariogramas empíricos pueden obtenerse usando la función variogram()

Variograma nube

v <- variogram(mg2040~ 1, 
                datos_sf_2, 
               cloud=T, 
               cutoff=700)
plot(v)

Variograma

semiv_emp <- variogram( mg2040~ 1, 
                       datos_sf_2, 
                       cutoff = 600)
head(semiv_emp)
##     np      dist    gamma dir.hor dir.ver   id
## 1  470  60.36988 26.33723       0       0 var1
## 2  656 108.13523 37.95655       0       0 var1
## 3  692 150.87729 43.01951       0       0 var1
## 4  575 181.80241 43.71304       0       0 var1
## 5  901 214.79382 44.33574       0       0 var1
## 6 1025 258.59180 46.49463       0       0 var1
plot(semiv_emp,
     main = "Contenido Magnesio (mmolc/dm3)",
     xlab = "Distancia",
     ylab = "Semivarianza")

En el caso anterior la fórmula utilizada (mg2040~1) asume que el proceso es estacionario.

A continuación se ajusta un modelo de semivariograma teórico sobre el semivariograma empírico

Esférico

mod_esf <- fit.variogram(
  semiv_emp,
  vgm(0.6, "Sph", 200, 0.2))
mod_esf
##   model    psill    range
## 1   Nug 10.41296   0.0000
## 2   Sph 36.06021 196.9126
plot(semiv_emp,
     mod_esf,
     main = "Contenido Magnesio (mmolc/dm3)",
     xlab = "Distancia",
     ylab = "Semivarianza")

El modelo que mejor ajusta será el de menor SCE. La función attr() devuelve atributos de un objeto y puede usarse para consultar la SCE del modelo ajustado.

attr(mod_esf, 'SSErr')
## [1] 0.4988919

Exponencial

mod_exp <- fit.variogram(
  semiv_emp, 
  vgm(0.6, "Exp", 200, 0.2))
mod_exp
##   model    psill    range
## 1   Nug  0.00000  0.00000
## 2   Exp 48.63929 75.68158
attr(mod_exp, 'SSErr')
## [1] 0.2319651
plot(semiv_emp,
     mod_exp,
     main = "Contenido Magnesio (mmolc/dm3)",
     xlab = "Distancia",
     ylab = "Semivarianza")

vgLine <- rbind(
  cbind(
    variogramLine(
      mod_exp, maxdist = max(semiv_emp$dist)),
    id = "Exponencial"),
  cbind(
    variogramLine(
      mod_esf, maxdist = max(semiv_emp$dist)),
    id = "Esférico")
  )


ggplot(semiv_emp, aes(x = dist, y = gamma, 
                      color = id)) +
  geom_line(data = vgLine) +
  geom_point() +
  labs(
    title = "Semivariograma experimental
    y teórico ajustado") +
  xlab("Distancia") +
  ylab("Semivarianza") +
  scale_color_discrete(
    name = "Modelo",
    labels = c("Esférico", 
               "Exponencial", 
               "Experimental"))

modelos <- fit.variogram(semiv_emp, 
                         vgm(c("Exp", "Sph" ,"Gau" ,"Mat")))
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
modelos
##   model    psill    range kappa
## 1   Nug  0.00000  0.00000   0.0
## 2   Mat 48.63929 75.68156   0.5
plot(semiv_emp,modelos)

attr(modelos, 'SSErr')
## [1] 0.2319651

Prediccion Espacial

La predicción espacial, es decir la predicción de valores de la variable en sitios del campo espacial donde no existen observaciones, usualmente se hace por el método kriging basándose en el semivariograma ajustado. Kriging proporciona el mejor estimador lineal insesgado del valor esperado para el sitio y un error de estimación conocido como varianza kriging. Esta varianza depende del modelo de semivariograma ajustado y de la ubicación en el espacio de los datos observados ya que son los datos observados en distintos sitios los que proveen información para aproximar el valor en el sitio sin dato.

limites <- read.csv("~/Regresion Avanzada/borders.csv", header=TRUE)
View(limites)
head(limites)
##   east north
## 1 5340  5800
## 2 5590  5690
## 3 5990  5690
## 4 5990  5100
## 5 5780  4800
## 6 5590  4800
grid <- pred_grid(limites, by = 10)
grid <- polygrid(grid, bor = limites)

names(grid) <- c("x", "y")
grid <- st_as_sf(grid, coords = c("x", "y"),
                 crs = 32721)
plot(grid)

En resumen, estos pasos permiten generar una grilla regular de puntos dentro de un polígono determinado y transformarla en un objeto que puede ser utilizado en análisis de estadística espacial.

La función krige() del paquete gstat se utiliza para realizar interpolación kriging y simulaciones condicionales con diferentes métodos de predicción. En este caso, se utiliza la interpolación por kriging ordinario con el modelo de semivariograma matern estimado con geoestadística clásica. La función toma como argumentos la fórmula para especificar la estacionariedad del proceso, la base de datos, el objeto para la predicción, y la información del modelo de semivariograma teórico ajustado. Los argumentos opcionales nmin y nmax permiten realizar la interpolación en un contexto local.

# Kriging Ordinario
kriging_o <-
  krige(
    mg2040 ~ 1,
    datos_sf_2,
    st_as_sf(grid),
    nmin = 7,
    nmax = 25,
    model = modelos
  ) 
## [using ordinary kriging]
#A continuación, se realiza la visualización de la predicción y su varianza.
predK_otm <- tm_shape(kriging_o) +
  tm_dots("var1.pred", style = "cont",
          title = "Predicción")
varK_otm <-  tm_shape(kriging_o) +
  tm_dots("var1.var", style = "cont",
          title = "Varianza")

tmap_arrange(predK_otm, varK_otm)

Validación cruzada

La validación cruzada k-fold es un proceso que se utiliza para evaluar la precisión de un modelo de interpolación. En esta secciónn realizaremos los siguientes pasos:

  1. Realizar la validación cruzada k-fold utilizando la función krige.cv() del paquete gstat.

  2. Utilizar los argumentos fórmula (mg2040~1, lo cual especifica que es un proceso estacionario), la base de datos (datos_sf_2) y el modelo de semivariograma teórico ajustado (modelos) en la función krige.cv().

  3. Definir el número de grupos (k) en los que se dividirá la base de datos utilizando el argumento nfold.

  4. Fijar la semilla mediante la función set.seed() para obtener repetibilidad en los resultados.

  5. Calcular estadísticos resumen como el error medio (ME), error cuadrático medio (MSE), media del cociente de la desviación cuadrática (mean squared deviation ratio, MSDR), raíz del error cuadrático medio (RMSE), la RMSE relativa a la media de los observados (RMSE_rel) y la correlación lineal entre los observados vs. predichos.

  6. Utilizar estos estadísticos resumen para evaluar el desempeño del modelo de interpolación kriging.

  7. Generar un gráfico de correlación entre los observados y predichos para visualizar la calidad del ajuste del modelo.

Realizar la validación cruzada, seleccion de k y escrbiri kriging

set.seed(17)
valcru <-
  krige.cv(
    mg2040 ~ 1,
    datos_sf_2,
    modelos,
    nfold = 10,
    nmin = 7,
    nmax = 25
  )

Calcular los estadisticos de resumen y predicciones

ME <- mean(valcru$residual)
MSE <- mean(valcru$residual ^ 2)
MSDR <- mean(valcru$zscore ^ 2)
RMSE <- sqrt(mean(valcru$residual ^ 2))
RMSE_rel <- 
  sqrt(mean(valcru$residual ^ 2)) / 
  mean(valcru$observed) * 100
r <- cor(valcru$observed, 
         valcru$observed - valcru$residual)

tabla <- data.frame(ME, MSE, RMSE, 
                    RMSE_rel, MSDR, r)
tabla
##           ME      MSE     RMSE RMSE_rel    MSDR         r
## 1 -0.2437268 27.86059 5.278313 20.65048 1.04758 0.6196965
plot(
  valcru$observed,
  valcru$observed - valcru$residual,
  xlab = "Observados",
  ylab = "Predichos"
)